NASA 

Technical 

Paper 

2098 

March 1983 


NASA 

TP 

2098 
c* 1 


Reduction Methods for lg|| 
Nonlinear Steady-State “■§ j 
Thermal Analysis 


Ahmed K. Noor, 
Chad D. Balch, 
and Macon A. Shibut 



NASA 




NASA 

Technical 

Paper 

2098 

1983 


TECH LIBRARY KAFB, NM 



0134 


as 


Reduction Methods for 
Nonlinear Steady-State 
Thermal Analysis 


Ahmed K. Noor, 

Chad D. Balch, 
and Macon A. Shibut 

The George Washington University 
Joint Institute for Advancement of Flight Sciences 
Langley Research Center 
Hampton, Virginia 


r\j/\s/\ 

National Aeronautics 
and Space Administration 


Scientific and Technical 
Information Branch 




CONTENTS 


SUMMARY 1 

INTRODUCTION 1 

MATHEMATICAL FORMULATION 2 

Governing Finite-Element Equations 2 

Solution of the Governing Equations 3 

BASIC IDEA OF REDUCTION METHODS 4 

Basis Reduction and the Reduced System of Equations 4 

Selection of Basis Vectors 5 

Comparison With Taylor Series Expansion 6 

Comments on the Selection of Basis Vectors 7 

COMPUTATIONAL ALGORITHM USED WITH REDUCTION METHODS 7 

Evaluation of the Basis Vectors and Generation of the Reduced System 

of Equations 8 

Sensing and Controlling the Error in the Reduced System of Equations .... 8 

NUMERICAL STUDIES 9 

Two-Dimensional Steady-State Conduction in a Square Plate 9 

Two-Dimensional Steady-State Conduction in a Cylinder With an 

Eccentric Hole 10 

Steady-State Analysis of a One-Dimensional Conducting- 

Convecting-Radiating Fin 11 

Steady-State Conduction and Radiation in a Segment of the Space 

Shuttle Orbiter Wing 11 

COMMENTS ON REDUCTION METHODS FOR NONLINEAR THERMAL ANALYSIS 12 

CONCLUDING REMARKS 13 

APPENDIX A - FORMULAS FOR COEFFICIENTS IN THE FINITE-ELEMENT EQUATIONS 

FOR INDIVIDUAL ELEMENTS 15 

APPENDIX B - EXPLICIT FORMS OF THE REDUCED SYSTEM OF EQUATIONS 17 

APPENDIX C - EVALUATION OF THE PATH DERIVATIVES 19 

APPENDIX D - TREATMENT OF PRESCRIBED BOUNDARY TEMPERATURES 26 

APPENDIX E - EXACT SOLUTION FOR TWO-DIMENSIONAL STEADY- STATE CONDUCTION 

IN A SQUARE PLATE WITH TEMPERATURE- DEPENDENT THERMAL CONDUCTIVITIES 28 

REFERENCES 29 

SYMBOLS 30 

FIGURES 33 


iii 



SUMMARY 


Hybrid analysis techniques based on the combined use of finite elements and the 
classical Bubnov-Galerkin approximation are presented for predicting nonlinear 
steady- state temperature distributions in structures and solids. In these hybrid 
techniques the modeling versatility of the finite-element method is preserved and a 
substantial reduction in the number of degrees of freedom is achieved by expressing 
the vector of nodal temperatures as a linear combination of a small number of global- 
temperature modes, or basis vectors. The Bubnov-Galerkin technique is then used to 
compute the coefficients of the linear combination (i.e., the amplitudes of the 
global- temperature modes) . 

The basis vectors chosen are the path derivatives commonly used in perturbation 
techniques, namely, the derivatives of the nodal-temperature vector with respect to a 
preselected control (or path) parameter or parameters. The vectors are generated by 
using the finite-element model of the initial discretization. Also, the performance 
of alternate sets of basis vectors is investigated. In the alternate sets, only a 
few path derivatives are generated, and they are augmented by a constant vector 
representing a uniform temperature rise (or drop) and by reciprocal vectors with non- 
zero components equal to the reciprocals of the nonzero components of the path deriv- 
atives. A problem- adaptive computational algorithm is presented for efficient evalua- 
tion of global approximation vectors and generation of the reduced system of equations 
and for monitoring the accuracy of the reduced system of equations. 

The potential of the proposed reduction methods for the solution of large-scale, 
nonlinear steady-state thermal problems is also discussed. The effectiveness of 
these methods is demonstrated by means of four numerical examples, including conduc- 
tion, convection, and radiation modes of heat transfer. 

This study shows that the use of the uni form- temperature mode and the path 
derivatives as global approximation vectors significantly increases the accuracy of 
the solutions obtained by reduction methods , thereby enhancing the effectiveness of 
these methods for the solution of large-scale, nonlinear thermal problems. 


INTRODUCTION 

Computational methods for nonlinear heat transfer have recently become the focus 
of intense research efforts because of the need for realistic modeling and accurate 
thermal analysis of large, complex hardware systems subject to harsh environments 
(e.g., reentry flight vehicle structures, large-area space structures, and nuclear 
reactor components) . Considerable progress has been made in the development of 
numerical discretization techniques. (See, for example, refs. 1 to 3.) Also, a 
number of versatile and powerful finite-element and finite-difference (lumped- 
parameter) programs have evolved for nonlinear thermal analysis. A survey of some of 
these programs is given in reference 4. In spite of these advances, the nonlinear 
thermal analyses of most large and complex hardware systems require excessive amounts 
of computer time even on present-day large computers and thus are very expensive. 

The large numbers of degrees of freedom required in the nonlinear thermal analy- 
sis of complex systems are often dictated by the topology of the system rather than 



by the expected complexity of the behavior. A similar situation has been observed in 
nonlinear structural analysis problems, and an effective reduction technique has been 
developed to significantly reduce the computational effort. (See refs. 5 and 6.) 

The application of the reduction method of reference 6 to transient nonlinear thermal 
problems is described in reference 7. The aim of the present study is to develop a 
reduction method and a computational algorithm for nonlinear steady-state thermal 
analysis of structures and solids. The proposed technique is similar to that pre- 
sented in references 5 and 6 for the nonlinear static analysis of structures and is a 
hybrid method which combines the modeling versatility of contemporary finite elements 
with the reduction in the total number of degrees of freedom provided by the classi- 
cal Bubnov-Galerkin technique. 

To sharpen the focus of the study, discussion is limited to nonlinear steady- 
state thermal problems with continuous temperature fields in the space domain. Con- 
vection, nonlinear conduction, and radiation modes of heat transfer are considered. 


MATHEMATICAL FORMULATION 
Governing Finite-Element Equations 

The solid region is discretized by using a single-field finite-element model with 
the fundamental unknowns consisting of the temperatures at the various nodes. The 
governing finite-element thermal-equilibrium equations can be cast in the following 
form: 


{ f (T) } = [K(T)]{T} - q{Q> = 0 


( 1 ) 


where {t} is the vector of nodal temperatures, {q} is the normalized thermal- 
load vector, [k(T)] is the n x n heat- transfer matrix (n is the total number of 
degrees of freedom) , and q is a normalizing parameter for the thermal-load vector. 
(A list of symbols used in this paper appears after the references.) The matrix 
[k (T) and the vector { q } can each be decomposed into three components as follows: 


[K(T)] = [K (k) ] + [K (r) ] + [ K < h >] (2) 

and 

(Q) = {Q (o) > + {Q (r) > + {Q (h) > (3) 

where [ K < k) ], [R (r) ], and [ K < h >] are the conduction, radiation, and convection 

matrices? and {Q^}, {Q^K and (Q ^ } are normalized applied heating, radia- 
tion, and convection thermal-load vectors. For convenience, the material conductivi- 
ties in the present study are assumed to vary quadratically with temperature, as 
follows : 



+ y T + Y 



) 


(4) 


2 



r 


where fe°g is the thermal-conductivity coefficient at a preselected reference 

temperature of T = 0 and y^ and y 2 are conductivity coefficients. In the 

numerical studies, only isotropic materials are considered, for which 

T + 

1 ± 

assumed to be independent of temperature and the radiation matrix varies cubically 
with temperature. 

The matrices [k^)] and [K (r) ] can be expressed in the following forms : 

(5) 

( 6 ) 


[K< k) ] = [K<°>] + qi [K (1) (T)] + q 2 [K (2) (T)] 

and 

[K< r >] = q 3 [K( r >] 


? T^ and k -^2 = 0. Also, the convection matrix is 


^11 ^22 ~ ^ i 1 + Y 


[K (o) ] is the linear conduction matrix (independent of T) ; [k^] and 

are normalized nonlinear conduction matrices whose terms are linear and 
quadratic in the nodal temperatures {t}; [k^] is a normalized radiation matrix; 

and q^, q 2 , and q^ are normalizing parameters. The expressions for [k^°^], 

[k (1) ], [k ( 2) ]/ [K (r) ]/ [K (h) ], {Q (r) >, and {Q (h) } for the individual 

elements are given in appendix A. 


where 

[k (2) ] 


Solution of the Governing Equations 

The solution of equation (1) is obtained by using an incremental-iterative tech- 
nique (e.g., a predictor-corrector continuation method). This is accomplished by 
embedding equation (1) in a single- or multiple-parameter family of equations of the 
form { f (T , q) } = 0 or {f(T, q-^ , q 2 , q 3 ) } = 0. Only the single-parameter case is 
discussed in this section. The multiple-parameter case is examined in subsequent 
sections (in conjunction with reduction methods) . The normalizing parameters q, 
q l' g 2' and q 3 are also referred to as control or path parameters. Henceforth, 
the terms normalizing parameter, control parameter, and path parameter are used 
interchangeably . 

In the single-parameter case, q^ = 3 2 = g 3 = For each value of q in some 

interval (e.g., 0 - q - q) , a solution {T(q)} exists which varies continuously 

with q and satisfies the conditions {t(0)} = {o} and {T(q)} = {t}. The solu- 
tion {t} corresponding to a particular value of the parameter q is used to calcu- 
late a suitable approximation (predictor) for {t} at a different value of q. This 
approximation is then chosen as an initial estimate for {t} in a corrective- 
iterative scheme such as the Newton- Raphs on technique. The iteration process is 
represented by the following equations for the j th iteration cycle: 


[j]<i>{AT}<i> = -{f(T,q)}(i> 


(7) 


3 



and 


{ T } ( i +1) _ { T }0) + {At} ^ ^ 


( 8 ) 


where [J] is the Jacobian matrix defined as 


[j] = [k] + 


8K I L t 
3t, T L 


(9) 


Note that the Jacobian matrix [j] is, in general, nonsymmetric . In the present 
study , in order to avoid the complications associated with solving nonsymmetric equa- 
tions, a quasi-Newton method is used. The method is based on using a symmetric 


approximation to the Jacobian matrix consisting of the sum of the matrix [k^°^] and 

^SK-t 

the contribution of the radiation matrix to 


"I L 


9t- 


The contributions of the 


nonlinear conduction matrices to 
neglected. 


9k 


IL 


3t, 


form a nonsymmetric matrix and are 


BASIC IDEA OF REDUCTION METHODS 
Basis Reduction and the Reduced System of Equations 

The essence of reduction methods is to reformulate the problem in terms of a 
few discrete variables {ij;}, which represent amplitudes of global-temperature modes 
(functions of the variables {t}) . This is accomplished by using the following 
trans f ormati on : 


{t} 


n, 1 




(r « n) (10) 


where [r] is an n x r transformation matrix in which the columns represent global 
approximation (or basis) vectors. Note that the number of generalized discrete 
variables r is assumed to be much smaller than the number of degrees of freedom of 
the initial discretization n. 

A Bubnov-Galerkin technique is then used to approximate equation (1) with a 
much smaller system of nonlinear algebraic equations in the new unknowns {ip } . The 
reduced equations have the following form: 


{f(T|>)} = [KWlW - q{Q} = 0 (11) 

where 

5 W] = [r] T [K(^)][r] (12) 


4 


and 


{Q} = [r ] T {Q} (13) 

The tilde (~) refers to the reduced system, the superscript T denotes transposi- 
tion, and [K(ip)] is obtained from [k(T)] by replacing {t} by its expression in 
terms of {if;} from equation (10) . 

The reduced-system heat- transfer matrix [k] and thermal-load vector {q} can 
be separated into components as follows (see eqs. (2), (3), (5), and (6)): 

[K(*J] = [K (0) ] + qi [K (1) ] + q 2 [K (2) ] + q 3 [K (r >] + [ic (h) ] (14) 

and 

{Q> = {Q (0) > + {Q (r) > + {Q (h) > (15) 

where [k^°^] is the linear conduction matrix; [k^] and [k^)] are the 

normalized nonlinear conduction matrices; [ K (h)] anc j [k^] are the convection 

and normalized radiation matrices; {Q^°^} is the normalized applied heating vector; 

and {q^} and {q^} are normalized radiation and convection thermal-load 
vectors. The explicit forms of the reduced-system arrays are given in appendix B. 

The solution of the reduced- system nonlinear equations (eq. (11) ) is obtained 
by using an incremental- iterative technique similar to the one described in the pre- 
ceding section for the full system of equations. However, since the number of equa- 
tions is small and their solution time is a small fraction of the total analysis 
time, no approximation needs to be made for the reduced-system Jacobian matrix [j] , 

and the nonsymme trie contributions of the nonlinear conduction matrices to 
are included in [j]. 


gi £ 


'Pc 


Selection of Basis Vectors 


The crux of reduction methods is the proper selection of the reduced basis 
vectors (columns of the matrix [r] in eq. (10)). An ad hoc or intuitive choice 
may not lead to a satisfactory approximation. A large number of numerical experi- 
ments with structural and solid mechanics problems have demonstrated that, in the 
case of a single path parameter (such as q) , the various-order derivatives of {t} 
with respect to q (path derivatives) provide an effective set of basis vectors. 

(See refs. 5 and 6.) In the present study, the effectiveness of using such basis 
vectors for thermal problems is investigated. The path derivatives are evaluated at 
q = 0 (q^ = q 2 = q^ = 1) and are given by 


[r] 





(16) 
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In the case of multiple path parameters (such as q^, q 2 , and q 3 ) , the basis 

vectors are chosen to be the various-order derivatives with respect to these path 
parameters evaluated at q^ = q 2 = <3 3 = 0 (q = 1 ) , that is, 


[r] 





2 





3 



(17) 


where {t} q is the linear solution (zeroth-order derivative) . The equations used in 
evaluating the basis vectors (path derivatives) are obtained by successive differen- 
tiations of the governing finite-element equations of the initial discretization in 
equation (1) . The explicit forms of these equations are given in appendix C for both 
the single-parameter and the multiple-parameter cases. The advantage of computing 
the path derivatives at zero values of the path parameters is that the left-hand-side 
matrix used in evaluating the path derivatives is independent of temperature. (See 
appendix C . ) 


The chosen set of basis vectors is the same as that commonly used in classical 
perturbation techniques. The basis vectors have the following properties: 


1. They are linearly independent and span the space of solutions in the neigh- 

borhood of the point of their generation, and therefore they fully charac- 
terize the nonlinear solution in that neighborhood. 

2. Their generation, using the finite-element model of the initial discretiza- 

tion, requires only one matrix factorization of the linear matrix 

[[k ( o) ] + [k (h) ]] (see appendix C) , and therefore it is computationally 
inexpensive . 

3. They provide a direct measure of the sensitivity of the thermal response to 

changes in the control (path) parameters. 


The first property is necessary for the convergence of the Bubnov-Galerkin 
approximation. The second property significantly enhances the efficiency of the 
reduction method and increases its effectiveness in solving large-scale nonlinear 
thermal problems. The implication of the third property is that by appropriate 
choice of the control parameters, sensitivity of the temperature distribution to 
changes in the thermal data of the medium (e.g., nonlinear conduction and radiation 
coefficients) can be obtained. Note that the sensitivity information is obtained at 
zero value (s) of the control parameter (s) . If the sensitivity is required at other 
values of the parameter (s) , approximate Taylor series expansion (s) may be used. 

It should be noted that the use of path derivatives as basis vectors in non- 
linear structural and solid mechanics problems provided highly accurate solutions 
within a large neighborhood of the point of evaluation of these derivatives 
(refs . 5 and 6 ) . 


Comparison With Taylor Series Expansion 

If the reduction method outlined in the preceding section is contrasted with 
the Taylor series expansion used in classical perturbation techniques, the following 
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can be noted. In both methods the vector of nodal temperatures {t} is approximated, 
over a range of values of the path parameter (s) , by a linear combination of the 
various-order derivatives of {t} with respect to the path parameter (s ) . However, 
the coefficients {tp} of the linear combination in the Taylor series expansion are 
fixed and are equal to 


(Aq) 2 (Aq) 

21 ' 3! ' 


for the single-parameter case and 


(Aq^ ) ^ 

1, Aq^, Aq 2 9 Aq 3 ' t Aq-^ Aq 2 , Aq^ Aq^, ... 

for the three-parameter case. By contrast, the coefficients {\p} in the reduction 
technique are left as free parameters and are determined by using the Bubnov-Galerkin 
technique. Therefore, the reduction method can be thought of as either of the 
following: 

1. A generalized Taylor series (or perturbation) approach with free parameters 

or coefficients 

2. A generalized Bubnov-Galerkin approach with the approximation (or basis) 

vectors generated by using a preselected finite-element model rather than 
chosen a priori 

Numerical experiments have shown that the use of the free parameters in the general- 
ized Taylor series technique leads to accurate solutions not only within the radius 
of convergence of the Taylor series but also well beyond it. (See refs. 5 and 6.) 


Comments on the Selection of Basis Vectors 

The computational effort required to generate the basis vectors can be reduced 
by generating only a few path derivatives (four or five) using equations (Cl) or (C2) 
and then augmenting these derivatives with a constant vector with equal components 
representing a uniform temperature mode, and/or reciprocal vectors with nonzero com- 
ponents equal to the reciprocals of the nonzero components of the vectors in equa- 
tions (Cl) or (C2) . The performance of the augmented set of basis vectors is dis- 
cussed in the section entitled "Numerical Studies." 


COMPUTATIONAL ALGORITHM USED WITH REDUCTION METHODS 

To realize the potential of reduction methods in large-scale, nonlinear thermal 
analysis, a problem- adaptive computational algorithm is needed which is both robust 
and efficient. The two key elements of the algorithm for nonlinear steady-state 
problems are the following: 

1. Efficient evaluation of the basis vectors and generation of the reduced sys- 

tem of equations 

2. Sensing and controlling the error in the reduced system of equations 
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These elements have been discussed in references 5 and 6 in connection with struc- 
tural and solid mechanics applications. The key features of the procedure are dis- 
cussed hereinafter. 


Evaluation of the Basis Vectors and Generation 
of the Reduced System of Equations 

The particular choice of the basis vectors to be the various-order derivatives 
of the nodal- temperature vector {t} with respect to the path parameter (s) permits 
the generation of all the vectors with only one factorization of the matrix 

[[K (0) ] + [K (h) ]] (See appendix C.) Therefore, the effort to generate the second 
and succeeding basis vectors reduces to that of evaluating the right-hand sides 
of equations (Cl) or (C2) for the single-parameter or the multiple-parameter cases. 
The expressions of the right-hand sides grow in complexity for higher order deriva- 
tives, and their computation involves contractions of multidimensional arrays with 
the basis vectors. (See appendix C.) To improve the computational efficiency, the 
contracted arrays which are common to more than one right-hand side are formed once 
and stored for subsequent use. 

Note that because of the recursive nature of the formulas for the path deriva- 
tives in both the single-parameter and the multiple-parameter cases, all the lower 
order derivatives must be evaluated before any subsequent derivatives can be com- 
puted. Once the entire set of path derivatives has been generated, it can be aug- 
mented with a constant vector, reciprocal vectors, or the vector {z} (eq. (D2) in 
appendix D) for the case of prescribed nonzero nodal temperatures to form the 
matrix [T]. 

The formation of the reduced arrays involves contraction of the full-system 
arrays with the basis vectors and appears to be the most time-consuming part of the 
solution process based on the reduction technique. The independent elements of the 
reduced arrays are generated once and stored for subsequent use, when the tempera- 
ture distribution is required for different value (s) of the path parameter (s) . 


Sensing and Controlling the Error in the 
Reduced System of Equations 

To check the accuracy of the solution obtained with the reduced system of equa- 
tions at any value of the path parameter, the approximate nodal temperatures are 
generated by using equation (10) with the vector {i/j } obtained by solving the 
reduced equations in equation (11) . Then the residual vector {r} of the original 
finite-element equations (eq. (1)) is computed as follows: 

{r} = [k(T)]{t} - q{Q> (18) 


In the present study, a weighted Euclidean norm of {r} is used as an error 
measure e, namely 


e = ^v/{R} T {R}/({Q} T {Q}) 

Si 


(19) 
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Numerical experiments have demonstrated that the error norm e is useful in 
assessing the overall accuracy of different reduction methods (based on different 
sets of basis vectors) . However, the magnitude of this error norm was found to be 
problem dependent, and therefore it is not possible to prescribe a value for the 
error tolerance which is suitable for different classes of nonlinear thermal problems. 
An alternate way of assessing the accuracy of the reduced solution is to use the 
vector of nodal temperatures {t} generated by the reduced system as a predictor in 
the full system, obtain a corrected estimate of {t} by performing a single iteration 
of the Newton-Raphson technique, and compare the predicted and corrected values 
of {t}. 


NUMERICAL STUDIES 

To test and evaluate the effectiveness of the proposed reduction techniques, 
several nonlinear thermal problems were solved. For each problem, solutions based 
on the full system of equations of the finite-element model were compared with other 
numerical approximations and exact solutions, whenever available; then the solutions 
obtained by the reduction method were compared with the full-system solutions and with 
other numerical approximations. The results of four typical problems are discussed 
herein. The four problems are the following: (a) two-dimensional steady-state con- 

duction in a square plate; (b) two-dimensional steady-state conduction in a cylinder 
with an eccentric hole; (c) steady-state analysis of a one-dimensional conducting- 
convecting-radiating fin; and (d) steady-state conduction and radiation in a segment 
of the Space Shuttle orbiter wing. In all the problems considered, the thermal con- 
ductivity is assumed to be temperature dependent. In the first two problems only the 
single-parameter reduction method is used, and in the last two problems both the 
single-parameter and the multiple-parameter methods are applied. 


Two-Dimensional Steady-State Conduction in a Square Plate 

The first problem considered is that of steady-state thermal conduction in a 
thin square plate with prescribed boundary temperatures. (See fig. 1.) The thermal 
conductivities vary quadratically with temperature. The material and geometric char- 
acteristics of the plate are given in figure 1. The problem was solved with a number 
of techniques, including an exact analytic technique based on the Kirchhoff trans- 
formation, finite elements, single-parameter reduction methods, and a perturbation 
technique based on Taylor series expansion of the temperature. The details of the 
exact analytic solution are given in appendix E. For the case of linear variation of 
conductivity coefficients with temperature, an exact analytic solution and a perturba- 
tion solution are presented in reference 8. 

Because of symmetry, only one-half of the plate was considered and was modeled 
by a grid of 5 x 10 elements with biquadratic Lagrangian interpolation functions for 
the temperature (a total of 190 temperature degrees of freedom (D.O.F.)). The 
finite-element solution was found to be in close agreement with the exact analytic 
solution. In both the single-parameter reduction method and the Taylor series 
expansion, the control parameter q was taken to be the amplitude of the nonzero 
prescribed boundary temperatures. Contour plots for the temperatures at q = 300 K, 
600 K, and 2400 K are shown in figure 1. As can be seen from figure 1, at higher 
values of q a boundary layer (with steep temperature gradients) is formed near the 
edges where a temperature of 0 K is prescribed. The basis vectors were evaluated at 
q = 0 K, and they were obtained by solving a linear set of finite-element equations. 

(See appendix C.) 


9 



Four sets of basis vectors are considered in the reduction method. The first 
two sets of basis vectors consist of the following: (a) the first eight derivatives 

with respect to the control parameter q; and (b) the first four derivatives with 
respect to q and four reciprocal vectors (with nonzero components equal to the 
reciprocals of the nonzero components of the path derivatives) . The last two sets of 
basis vectors are the same as the first two, except for the addition of a constant 
vector. Normalized contour plots for the first six path derivatives are shown in 
figure 2. 

The accuracy of the solutions obtained by using the Taylor series expansion and 
the reduction method is indicated in figures 3 and 4. Also, figure 5 shows the vari- 
ations of the error norm e (eq. (19)) and the root-mean- square error as the number 
of reduced-system basis vectors increases. An examination of the results shown in 
figures 3 to 5 reveals the following: 

1. The Taylor series solution is considerably in error, even for the lowest value 

of q (q = 300 K) . (See fig. 3.) The radius of convergence of the Taylor 
series was found to be in the neighborhood of q = 240 K. By contrast, the 
solutions obtained by using the single-parameter reduction method with eight 
basis vectors were generally accurate, except near the boundary zone. (See 
fig. 4.) 

2. The accuracy of the solutions obtained by using four derivatives and four 

reciprocal vectors is comparable to that obtained by using the eight deriva- 
tives. (See fig. 4.) However, the computational effort involved in gen- 
erating the eight derivatives is considerably more than that required for 
generating the four derivatives and the reciprocal vectors . 

3. The accuracy of the solutions obtained by using the reduction method was con- 

siderably improved when the eight basis vectors were augmented by a constant 
vector. This is true for both sets of basis vectors. (See fig. 4.) 

4. The variation of the error norm e is similar to that of the root-mean- square 

error (see fig. 5) , and the solutions obtained with the reduction method 
converge to the full-system solution as the number of basis vectors 
increases . 


Two-Dimensional Steady-State Conduction in a Cylinder 
With an Eccentric Hole 

The second problem is steady-state conduction in a cylinder with an eccentric 
hole. The temperatures of the inner and outer surfaces of the cylinder are prescribed 
to be 1000 K and 0 K. The finite-element model along with the material and geometric 
characteristics of the cylinder is shown in figure 6. A previous finite-element 
solution is presented in reference 9 for linear variation of the conductivity coeffi- 
cients with temperature. Because of symmetry, only one-half of the cross section was 
modeled by a grid of 48 isoparametric, 9-noded elements (a total of 195 temperature 
degrees of freedom) . Contour plots for the temperatures are shown in figure 6 for 
the two cases of y 2 = 0 and Y 2 = which correspond to linear and quadratic 

variation of thermal conductivity with temperature. A boundary layer (with steep 
temperature gradients) exists near the outer surface when y 2 = 0.005, as shown in 
figure 6. The control parameter q was taken to be the amplitude of the prescribed 
boundary temperature. The basis vectors were generated at q = 0. 
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Four sets of basis vectors are again considered, namely, eight derivatives with 
respect to the control parameter q, eight derivatives and a constant vector, four 
derivatives with respect to q and their reciprocal vectors, and four derivatives, 
their reciprocal vectors, and a constant vector. The accuracy of the solutions 
obtained by using each of the sets of basis vectors is indicated in figure 7 for the 
two cases of Y 2 = 0 and y 2 = 0.005. As in the previous problem, the addition of 
the constant vector considerably improves the accuracy of the solutions, particularly 
near the outer surface for the case of y 2 = 0.005. (Note that the error in the solu- 
tion obtained without the use of the constant vector is represented by the distance 
between the cross and circle symbols in fig. 7.) Also, the accuracy of the solutions 
obtained by using four derivatives and four reciprocal vectors is comparable to that 
obtained by using the eight derivatives. 


Steady-State Analysis of a One-Dimensional 
Conducting-Convecting-Radiating Fin 

The third problem considered is steady-state analysis of a one-dimensional 
conducting-convecting- radiating fin. (See fig. 8 .) The temperature at one end is 
prescribed to be 1000 K. The fin was modeled by using a variable grid of 15 one- 
dimensional finite elements having 3 nodes and quadratic Lagrangian interpolation 
functions for the temperature (a total of 30 temperature degrees of freedom) . The 
thermal conductivity is assumed to vary linearly with temperature, and therefore the 
nonlinear conduction matrix [k ^ (T) ] is zero. 

Solutions were obtained by using both the single-parameter and the two-parameter 
reduction methods. In the single-parameter method, the control parameter q is the 
magnitude of both the prescribed temperature at x^_ = 0 and the convection thermal- 
load vector. The first six derivatives with respect to q were generated at q = 0. 
In the two-parameter method, the control parameters are chosen to be q 1 and q^ . 

Six path derivatives were generated at q^ = q^ = 0. (These include the linear solu- 
tion and all the first and second derivatives with respect to q^ and 33 *) The 
accuracy of the solutions obtained by using the reduction methods is indicated in 
figure 9. The solutions obtained by using both the single-parameter and the two- 
parameter reduction methods are in close agreement with the full-system solution 
except near the edge where x^ = 1.0. The maximum errors in the temperature at 
xj_ = 1.0 obtained with the single-parameter and the two-parameter reduction methods 
are 6.9 and 5.8 percent. In the presence of convection, it was necessary to add the 
constant vector to the two sets of basis vectors to obtain accurate solutions. 


Steady-State Conduction and Radiation in a Segment 
of the Space Shuttle Orbiter Wing 

The fourth problem is steady- state conduction and radiation in a segment of the 
Space Shuttle orbiter wing (fig. 10) . The problem was selected to assess the effec- 
tiveness of reduction methods for complex aerospace applications. 

The model consists of a two-dimensional cross section of a single bay of the 
Shuttle orbiter wing including reusable surface insulation (RSI) and RSI coating, 
room temperature- vulcanizing (RTV) rubber, a strain isolation pad (SIP) , and part 
of the aluminum wing box structure. This model is an adaptation of a finite-element 
model developed by the authors of reference 10. The applied heating on the structure 
consists of a distributed heat source on the outer surface Q(xj_), which varies 
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quadratically with x The temperature of the aluminum inner surface is prescribed 
as 394.4 K and the outer surface is a radiating boundary. The material and geometric 
characteristics of the structure and the finite-element model used are shown in fig- 
ure 10. Note that the scale of the model is magnified in the x 2 direction for 
visual clarity. Rectangular nine-noded Lagrangian elements with biquadratic interpo- 
lation functions for the temperature are used in modeling the RSI , RTV, and SIP. The 
aluminum structure, RSI coating, and radiating boundary are modeled by using one- 
dimensional three-noded elements with quadratic temperature interpolation functions. 
The finite-element model has a total of 64 nine-noded rectangular elements, 16 three- 
noded conduction elements, and 8 three-noded radiation elements (a total of 272 tem- 
perature degrees of freedom) . 

Solutions were obtained by using both the single -parameter and the three- 
parameter reduction methods. In the single-parameter method, the control parameter q 
was the magnitude of both the prescribed temperature and the applied heating. A total 
of eight path derivatives were generated at q = 0. In the three-parameter case, 

10 path derivatives were generated at q^ = q 2 = q 2 = 0. 

The following two sets of basis vectors were used: the first eight derivatives 

with respect to the control parameter, and the linear solution and all the first- and 
second-order derivatives with respect to the control parameters q^, q 2 , and q^. A 
constant vector was added to each of the two sets. The accuracy of the solutions 
obtained by using each of the two sets of basis vectors is indicated in figure 11, in 
which plots of the temperature distribution along lines parallel to the x-j_ and x 2 
axes are shown. The variations of the error norm e and the root-mean-square error 
with the number of basis vectors are shown in figure 12. The addition of a constant 
vector to the path derivatives was necessary to obtain acceptable accuracy. (Results 
obtained without the use of the constant vector are not shown.) 

An examination of the results shown in figures 11 and 12 reveals the following: 

1. Both the single-parameter and the three-parameter reduction methods give 
accurate results for the high temperatures at the outer surface. However, the solu- 
tions obtained by using the single-parameter method are more accurate than those 
obtained by using the three-parameter method in the low- temperature region (near the 
aluminum structure) . 

2. The variation of the error norm e is similar to that of the root-mean-square 
error, and the solutions obtained with both the single-parameter and the three- 
parameter reduction methods converge to the full-system solution as the number of 
basis vectors increases . 
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COMMENTS ON REDUCTION METHODS FOR NONLINEAR THERMAL ANALYSIS 

The proposed reduction methods appear to have high potential for solution of 
large-scale, nonlinear thermal problems. The numerical studies conducted clearly 
demonstrate the accuracy and effectiveness of the techniques. In particular, the 
following two points are worth mentioning: 

1. The small size of the reduced system of equations makes it feasible to con- 
duct parametric studies by repeatedly solving the equations for varying values of 
the control parameter (s) . 
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2 . For the problems considered, the reduction in the analysis time achieved by 
using the proposed technique in general was not proportional to the reduction in the 
number of degrees of freedom of the initial discretization. This was due to the rela- 
tively high cost of evaluating the basis vectors and of generating the reduced equa- 
tions (as compared with the cost of solving the reduced nonlinear algebraic equa- 
tions) . For any given problem, the saving in the analysis time resulting from the 
use of reduction methods is a function of the hardware and of the system software as 
well as the efficiency of the program used in the analysis, and therefore it is diffi- 
cult to quantify. In the numerical examples considered, overall reductions in the 
CPU times are estimated to be factors of 2 to 3. However, it should be noted that the 
program used for the computation was not optimized, and additional information in the 
form of sensitivity of the thermal response was provided by the basis vectors used in 
the proposed technique. More dramatic savings in CPU times are expected in large- 
scale problems with thousands of degrees of freedom if an efficient computer imple- 
mentation of the proposed technique is made. 


CONCLUDING REMARKS 

Reduction methods and a problem-adaptive computational algorithm were presented 
for predicting nonlinear steady-state temperature distributions in structures and 
solids. The computational algorithm can be conveniently divided into the following 
two distinct stages: discretization by using the finite-element method, and reduction 

in the number of degrees of freedom of the initial discretization by expressing the 
vector of unknown nodal temperatures as a linear combination of global-temperature 
modes, or basis vectors. A Bubnov-Galerkin technique is used to compute the coeffi- 
cients of the linear combination (amplitudes of the global- temperature modes) . The 
basis vectors are chosen to be those commonly used in the perturbation technique, 
namely, the derivatives of the nodal- temperature vector with respect to preselected 
control (or path) parameters. The vectors are generated by using the finite-element 
model of the initial discretization. 

Four numerical examples demonstrate the effectiveness of reduction methods for 
the solution of nonlinear steady-state thermal problems. The four problems are two- 
dimensional steady-state conduction in a square plate, two-dimensional steady-state 
conduction in a cylinder with an eccentric hole, steady-state analysis of a one- 
dimensional conducting-convecting-radiating fin, and steady-state conduction and 
radiation in a segment of the Space Shuttle orbiter wing. In all the problems con- 
sidered, the thermal conductivity is assumed to be temperature dependent. 

The results of the study suggest several conclusions relative to the selection 
of basis vectors and to the effectiveness of using reduction methods in nonlinear 
steady-state thermal problems. These conclusions are as follows: 

1. The proposed reduction methods are hybrid techniques which combine the major 
advantages of the finite-element numerical discretization technique, the classical 
Bubnov-Galerkin technique, and the perturbation method. These advantages are model- 
ing versatility, reduction in total number of degrees of freedom, and simplicity of 
assessing the sensitivity of the thermal response to variations in the control (or 
path) parameters. Moreover, the proposed methods greatly alleviate the following 
major drawbacks of the aforementioned three techniques : 

a. Excessive amounts of computer time required for the nonlinear finite-element 
analysis of large, complex systems 
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b. Difficulty of selecting global approximation functions for the classical 

Bubnov-Galerkin technique 

c. Small radius of convergence of the Taylor series expansions used in the 

classical perturbation techniques 

2. The use of path derivatives as basis vectors in nonlinear steady-state 
thermal problems leads to accurate solutions with a small number of basis vectors. 
Therefore, the time required to solve the reduced equations is relatively small, and 
the total analysis time to a first approximation equals the time required to evaluate 
the basis vectors and generate the reduced equations. 

3. The range of applicability of the reduction method can be greatly extended by 
augmenting the path derivatives of the nodal temperatures with a constant vector 
representing a uniform temperature rise or drop. In the problems considered, the 
addition of a constant vector was crucial for obtaining accurate solutions. 

4. The computational effort required for generating the basis vectors can be 
reduced by generating only a few path derivatives (4 or 5) and then augmenting these 
vectors with "reciprocal vectors," in which nonzero elements are equal to the recip- 
rocals of the nonzero values of the path derivatives. The accuracy obtained by using 
this set of basis vectors is comparable to that obtained by using an equal number of 
path derivatives. 

5. The accuracy of the multiple-parameter reduction method is comparable to that 
of the single-parameter method. However, the former provides more information about 
the sensitivity of the thermal response to variations in the material properties of 
the structure. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
December 6, 1982 
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APPENDIX A 


FORMULAS FOR COEFFICIENTS IN THE FINITE- ELEMENT 
EQUATIONS FOR INDIVIDUAL ELEMENTS 

The finite-element equations for individual elements can be written in the fol- 
lowing form: 

[ K Ij’ + '3l K MK T K + VljLVl + ’3 K S™ T K T L T N + K S’] T J - + ef’ + fif’) 

(Al) 

where , K^P., and are full-system linear and normalized nonlinear con- 

Iu IJK IJKL 

duction arrays; and are full-system radiation and convection arrays; 

Q^°^ are normalized applied heating thermal-load coefficients for the full system; 
and Q-[ r ^ and are normalized radiation and convection thermal-load coeffi- 

cients for the full system. The range of the uppercase subscripts (Latin index) is 
1 to m, where m is the number of nodes in the element, and a repeated uppercase 
subscript denotes summation over the full range. 


The coefficients in equation (Al) are defined as follows: 
.(o) 


K. 


= / k° Q 3 M t 3„N t d Q 

ij J a3 a i 3 J 


K 


(1) = 
IJK 


/ 


SI 


(e) 


Va3 d a N i Vj N k dQ 


.(2) 


= J y fe° Q 3 W 

1 O A/fi nj J 


K iA - J , Y 2 K aB Vi Vj n A dn 

n (e) 


/ 


(r) 

'IJKLN " (e) I J K L N 


= J aeNJlJl.JN_N„ dc 

I 

c 


.(h) 


= f hW w 


K = J hN N dc 
IJ (e) I J 

c h 
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. (o) 


f QW t an 


,(r) 


= f aHW 
(e) 1 


dc 


. (h) 


L 


hT W dc 
h I 


In the above equations , Mj are shape (interpolation) functions, 3 a = 3/3x a , 

( e ) ( 0 ) 

3g E 3/3xg, Q is the element domain, c^ is the portion of the element domain 

(e) 

where radiant heating (or cooling) is prescribed, c^ is the portion of the element 


domain where convective heating (or cooling) is prescribed, Q is the source heat 
rate, 0 is the Stefan-Boltzmann constant, £ is the emissivity, a is the absorp- 
tivity, H is the incident radiant heat flux, h is the convection coefficient, and 
T h is the known convective-exchange temperature. The range of the indices a, 3 is 
1,2, and a repeated index denotes summation over the full range. Note that the K 

arrays possess the following indicial symmetries: and K_^ are symmetric 

(i) 

with respect to interchange of I and J; K is symmetric with respect to mter- 

(2) 

change of I and J; K is symmetric with respect to interchange of I and J 

IJKL 

( 2r) 

or K and L; and K is completely symmetric with respect to interchange of 

I J KXiN 

any two indices. 


Elements with identical geometry and identical thermal properties have identical 
K arrays. Therefore, the K arrays are generated only for different elements 
(i.e., for elements with different geometry and/or different thermal properties). 
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EXPLICIT FORMS OF THE REDUCED SYSTEM OF EQUATIONS 

The reduced system of equations (eq. (11) ) can be written in index notation as 
follows : 

[SJ? + + + V&JtnWn ♦ ' -»(s‘ o) + 5{ r> + fif’) 

(Bl) 

The structure of these equations is similar to that of the finite-element equations 
for individual elements given in appendix A. The arrays in equation (Bl) can be 
expressed in terms of the full-system arrays defined in appendix A and the matrix of 
basis vectors F^ as follows: 


K 


(o) 

ij 


Z K (o) r .r . 

IJ II Jj 


Elements 


Jj 


K 


(i) _ K VJ_; r r r 

ijk J IJK Ii Jj Kk 


.(D- 


Elements 


t ^(2) __ \ S jr ( 2) -p -p p -n 

ijk£ Z-J IJKL Ii Jj Kk Li 

Elements 


= \ s r r r r r 

ijk£n Z-J IJKLN Ii Jj Kk l£ Nn 


(r) 


Elements 


K 


. (h) = K:“'r,.r 

1] 


Elements 


(h). 

IJ * Ii* J j 


~,o, 


Elements 


Q (o) r . 

V I Ii 


(r) 


Elements 


Q (r) r . 

V I Ii 


Qf>= 2 8j hl Tj, 

Elements 
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where the range of uppercase subscripts (Latin index) is 1 to m (number of nodes in 
the element) , the range of lowercase subscripts (Latin index) is 1 to r (number of 
reduced-system degrees of freedom) , and a repeated index in the same term denotes 
summation over the full range. 


Note that the indicial symmetries exhibited by the full-system arrays are 


retained in the reduced system as follows: 


:( 1 ) 


K (0) 
1 J 


and K. (h) 
ID 


K ijk s ^ inmetr ^ c respect to interchange of i 


and 


are symmetric arrays; 

• ?( 2 ) • _ . 
j; K _. 41 , o 1S symmetric 

with respect to interchange of i and j or k 
symmetric with respect to the interchange of any two indices. 


•V (r*) 

and Z; and K. n 
l jk£n 


is completely 
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EVALUATION OF THE PATH DERIVATIVES 

The path derivatives in equations (16) and (17) are obtained by successive dif- 
ferentiation of the governing finite-element equations (eq. (1) ) with respect to the 
path parameter q in the single-parameter reduction method and with respect to the 
parameters q^, q^/ and q^ in the three-parameter reduction method. The explicit 

forms of the recursion relations for individual elements are given in this appendix. 


Single-Parameter Reduction Method 

The recursion relations for the path derivatives can be written in the following 
compact form: 


[*£’ * ‘S’] 


3 a T- 




= R 


(a) 


(Cl) 


q=0 


where I and J range from 1 to m (number of nodes in the element) , a ranges 
from 1 to r (number of reduced-system degrees of freedom) , and the explicit forms 

of the components of the right-hand side R.f a ^ are given in table Cl. 


TABLE Cl.- EXPLICIT FORMS OF COMPONENTS OF R. 


(a) 


Explicit form of R 


FOR SINGLE-PARAMETER METHOD 
(a) 


Q (0) + Q (r) + Q (h) 
^1 *1 ^I 


-2K 


(1) 9T J 9t k 
IJK 3q 3q 


- 3K. 


(l) / 3 Tj 9 2 t k 3 2 Tj 3t k ' 


-6K 


i JK \3q 3q 2 3q 2 3q 

(2) 3t J 3t K 3t L 
IJKL 3q 3q 3q 


-K. 


(1) / 3 Tj 3 3 T k 3 2 Tj 3 2 t 


IJK 


3q 3 + 6 3q 2 3q 2 


K x . 3 3 Tj 9Tk 
+ 4 


3q 


3 3q 


-K 


( 2 ) 

'IJK L 


,24 


3 Tj 9t k 32t l 
3q 3q 9q 2 


+ 12 


3 2 Tj 9t k 9t l 


3q 


2 3q 3q 


-24K 


(r) 3Tj 3 T k 3t l 3T N 
IJKLN 3q 3q 3q 3q 
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TABLE Cl • ~ Continued 
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TABLE Cl.- Concluded 


Explicit form of R. 


(a) 


, K u) i 8 ps + 28 

IJK\ 9q 


9 6 t t 9 2 tv 9 5 t t 9 3 t„. 9 4 t t 9 4 t k 9 3 Tj 9 5 t k 

J K + 56 — J- — ~ + 70 r~ A +56 J — 


3q' 


9q 6 9q 2 


9q 5 9q 3 


+ 28 


8 2 Tj 9 6 t k 
9q 2 9q 6 


+ 8 


9Tj 9 7 T k \ 


9q 




-K 


( 2 ) 

I JKL 


56 


9 6 Tj 9t k 9t l . _ 95t j 9t k 92t l , 

+ 336 = — -r— z — + 


9 4 t 


9q 


6 9q 9q 


9q' 


9q 


9q 4 9q 4 


T , 9t k 9 3 t t 
J 560 K L 


3q' 


9q^ 


9q 9q 3 


d 2 i v 9 2 t t 

+ 420 — ^ + 


9q 9q^ 


9 3 T 


, , 9t„ 9 4 T t 9 2 Tt, 9 3 T t 

J 560 _* A + H20 K L 


9q- 


9q 9q 4 


9q^ 9q~ 


9 2 T- 


9t„ 9 5 T. 


3q 2 

t 36 31 

3q 

9t t 

/ 9t k 

112 TT 

9 6 t : 


j d 2 T„ 9 4 t t 9 3 t k 9 3 t t 

L + 840 — ^ ^ + 560 — ^ 

9q 3 9q 3 


9q z 9q 4 


9q 


9q 


9q fc 


, d 2 T v 9 5 t t 9 3 t k 9 4 t t 

^ + 336 — ^ ^ + 560 — ^ 

9q 3 9q 4 


9q 2 9q 5 


(r) 


( 9t k 9t l 9 5 t n 9 Tj 9t k 9 2 t l 9 4 t n 

' K IJKLN\^ 1344 Sq 9q 3q a q 5 + 3q 3q g q 2 3q 4 

3Tj 9 2 t k 9 2 t l 9 3 t n 9Tj 9t k 9 3 t l 9 3 t n 

+ 6720 


+ 20 160 


9< 2 9q 2 9q 2 


3q J 


9q 9q 


9q 3 9q 3 


+ 2520 


9 2 Tj 9 2 T K 9 2 t l 9 2 t^ 


9q z 9q" 


9q 2 9q 2 


9q 9q~ 


Three-Parameter Reduction Method 

The recursion relations for the individual elements have the following form: 


[*£' * <1 


,p+S+t n 


9q ? 3q S 2 3q^ 


= R. 


(a) 


(C2) 


q l =q 2 =q 3 =0 


where I and J range from 1 to m. The total number of p+s+t combinations 
is r, and the explicit forms of the components of the right-hand side are 

given in table C2. 
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,(a) 


TABLE C2 . - EXPLICIT FORMS OF COMPONENTS OF FOR THREE- PARAMETER METHOD 


Order 

of 

equation 


Zeroth 


First 


Second 


Explicit form of R, 


(a) 


Q (o) + Q ( r) + 

*1 *1 *1 


- K wi T J T K 


-K^^TTT 
XiJisJj J is. J_i 
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— K x 1 T T T T 
IJKLN J K L N 


- 2 K a> (- ^ 
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+ T 


IJKy J 3q^ K 
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IJKLl j K 3q 2 K L 3q„ 
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3 t t 
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TABLE C2.~ Continued 


Explicit form of 



. K d) L + 

IJKy J 3q 3 K dq_2 


(2) ( 9T T 9T 7 

1JKL y 1 J 1 K 3q 1 1 K i L 3q 1 


(1 ) / 3^t k 3t t 3t k 9 t t\ 

- K -( 3 T ^ + 6 ^r^ + 3T ^j 
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Explicit form of R 
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(r) / 92t n 3T L 9t n 

' K IJKLN ( 8T J T K T L 3q 1 3q 3 + 24T J T K 3 q;L 3q 3 


3T 3t 3 t k 

| 9 iL — + 2T — + 2 

K IJk 1 8q i 9q 2 J 9q i 9q 2 


e>\ / 3 t t 3t k 3t l 

K IJKL^“ T J T K ^2 J 9q x 9q 1 


3t 

3T„ 

3 2 t t \ 

J 

^ , om 


9q 2 

3q 2 k 

3q x 3q 2 j 


3T J 3t l 

3 2 t 

rrt m 


K 3q x 3q K L 2 


. K W/ T UljS+j— — + T ^ 
K IJ K | J 9q 2 3q 2 3q 2 K }q 2 
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In tables Cl and C2 a repeated uppercase Latin index denotes summation over the 
range 1 to m. The evaluation of the right-hand- side vectors listed in tables Cl 
and C2 involves contraction of K arrays with basis vectors. Many of the contracted 
arrays reappear in subsequent right-hand sides. For computational efficiency, the 
contracted arrays which are common to more than one right-hand side are generated 
once and stored for subsequent use. Note that the coefficient matrix on the left- 
hand sides of equations (Cl) and (C2) , which must be factored, is the same for each 
of the path derivatives. Hence, after assembly, the global matrix is factored only 
once regardless of the number of path derivatives generated. 
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TREATMENT OF PRESCRIBED BOUNDARY TEMPERATURES 

In the case of prescribed nonzero values for the boundary temperatures , it is 
convenient to partition the vector of nodal temperatures {t} as follows: 



where {T f } and {t } are the vectors of free and prescribed nonzero nodal tempera- 
tures. The prescribed zero nodal temperatures and their associated equations are 
eliminated from equation (1) . The prescribed nonzero nodal temperatures are assumed 
to be proportional to the parameter q, that is. 


{T p } = q{z} 


(D2) 


Governing Finite-Element Equations 

Equation (1) can be conveniently partitioned into two sets of matrix equations 
as follows: 


|f f (T f ,T p ) 


lyvy 


K £f (T) K fp (T) 


K pf (T) K pp m 



r ~\ 

% 




(D3) 


In the absence of heat sources and radiation and convection thermal loads, 

{ Q £ } = 0 and {Qp} equals the equivalent thermal loads associated with the pre- 
scribed temperatures {T }. The first set of equations in equation (D3) is used to 

ir 

determine {T f } and the second set can then be used to evaluate the vector (Qp}« 

Basis Vectors and the Reduced System of Equations 

in equation (10) can be conveniently partitioned as follows: 


(D4) 


The matrix [r] 


[r] = 
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where 
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for the single-parameter case, and 
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0 

V J 

... 


for the three-parameter case. 


(D5) 


(D6) 


The corresponding reduced-system equations are given by equation (11) with 
given by 




{Q> = 



(D7) 


Note that when {Qf} = 0, {Q} has only one nonzero component, namely . Equa- 
tion (11) is solved for the reduced unknowns {if;} subject to the condition ip 1 = q. 

Whenever desired, the basis vectors in equations (D5) and (D6) are augmented by 
the constant and reciprocal vectors defined in the paper. 
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EXACT SOLUTION FOR TWO-DIMENSIONAL STEADY-STATE CONDUCTION IN A 
SQUARE PLATE WITH TEMPERATURE-DEPENDENT THERMAL CONDUCTIVITIES 

The exact solution for the two-dimensional steady-state conduction in a square 
isotropic plate with temperature-dependent thermal conductivities can be obtained by 
using the Kirchhoff transformation (see page 11 of ref. 11) , which involves the fol- 
lowing substitution : 


U(T) 



fe(T) dT 


(El) 


For quadratic variation of conductivity with temperature of the form 
fe(T) = fe°^l + y^T + y 2 ^ T 2 )/ equation (El) becomes 

U(T) = T (l + | Y 1 T + j Y 2 T 2 ) 


(E2) 


The governing differential equation for heat conduction in the plate is linear 
in the transformed variable U(T). The expression of U(T) which satisfies the 
boundary conditions shown in figure 1 can be written in the following form: 


U 


oo 

2 : 


mTTx^ mTTx2 

A sin sinh * 

m l L 


(E3) 


where L is the side length of the plate and 


^m 


L sinh(mTr) 


f 


mTfXn 


sm 


7TX-, 


q sin 


hi(’ !l, T) + 7^(' 


q sm 


ttx-l \3 


dx-. 


(E4) 


Once the coefficients A m have been evaluated and U is determined from equa- 
tion (E3 ) , the value of the temperature T at each point can be obtained by solving 
the cubic equation (eq. (E 2) ) . 

In the present study, 10 terms were included in the series expansion of U 
(eq. (E3) ) , with an estimated maximum truncation error of 1 x 10 “ ^ for the tempera- 
ture T. The computerized symbolic and algebraic manipulation system MACSYMA 
(ref. 12) was used in the computation. 
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SYJyiBOLS 


A cross-sectional area of fin 

a absorptivity 

c circumference of fin 

{ e ) ( 0 ) 

c portions of the element domain where radiant and convective heating 

(or cooling) are prescribed (see appendix A) 

e error norm defined in equation (19) 

{f (T) } , { f (ip) } vectors defined in equations (1) and (11) 


H incident radiant heat flux 

h convection coefficient 

[j] Jacobian matrix defined in equation (9) 

[j] Jacobian matrix of the reduced equations 


[k(t)],[kojo] 

[K< k q,[K (r) ], 
[K<r)] [K<h>] 

[k (o) ],[k ( o) ] 

[k (1) ],[k (2) ], 

[k(D],[k( 2 )] 

[K< r) ] [K (h >] 


K (o) (1) (2) 

IJ ' IJK' IJKL 


K (r) k (h) 
IJKLN' IJ 

~(o) ~(1) ~ (2) 
ij ' ijk' ijk£ 


~(r) ~(h) 

ijk£n' ij 

k° 

k Q ,k° Q 
aB aB 


heat- transfer matrices of the full and reduced systems (see 
eqs. (1) and (12)) 

conduction, radiation, normalized radiation, and convection matrices 
of the full system (see eqs. (2) and (6)) 

linear conduction matrices (independent of T) of the full and 
reduced systems (see eqs. (5) and (14)) 

normalized nonlinear conduction matrices of the full and reduced 
systems (see eqs. (5) and (14)) 

normalized radiation and convection matrices of the reduced system 
(see eq. (14) ) 

full-system linear and normalized nonlinear conduction arrays 
defined in appendix A 

full-system radiation and convection arrays defined in appendix A 

reduced-system linear and normalized nonlinear conduction arrays 
defined in appendix B 

reduced-system radiation and convection arrays 

thermal-conductivity coefficient of an isotropic material 

thermal-conductivity coefficients in equation (4) 


30 



m 


number of nodes in an element 


n 


total number of temperature degrees of freedom in the analysis model 


bl shape (interpolation) function 

{Q},{Q} normalized thermal-load vectors for the full and reduced systems 


{QfMQp} 


thermal- load vectors associated with free and prescribed temperatures 
(see eq. (D3)) 


{Q (o) },{Q (r) },{Q (h) } 

{Q (o) },{Q (r) },{Q (h) } 


2 i 0) '2^ r) '2l h) 


Q. (o) ,Q. (r) ,Q. (h) 
1 1 


normalized applied heating, radiation, and convection thermal- 
load vectors for the full system (see eq. (3)) 

normalized applied heating, radiation, and convection thermal- 
load vectors for the reduced system (see eq. (15)) 

normalized applied heating, radiation, and convection thermal- 
load coefficients for the full system (defined in 
appendix A) 

normalized applied heating, radiation, and convection thermal- 
load coefficients for the reduced system (defined in 
appendix B) 


q,q^,q 2 ,q^ normalizing (or control) parameters 

{R} residual vector defined in equation (18) 


number of basis vectors (reduced- system degrees of freedom) 
temperature 


{t},{t £ },{ T p } 


vectors of free and prescribed nodal temperatures (see 
(eq. (Dl)) 


convective-exchange temperature 
component of nodal- temperature vector 
Cartesian coordinates 

vector of normalized prescribed nonzero nodal temperatures (see eq. (D2) ) 
conductivity coefficients (see eq. (4) ) 

[r]f[T f ],[r ] matrices of basis vectors defined in equations (10) and (D4) 

±r 

£ emissivity 

{ip} vector of undetermined coefficients of the reduced system 

ip. undetermined coefficients of the reduced system (amplitudes of global- 

temperature modes) 


x l' x 2 

{z} 

Yl'Y 2 
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a Stefan- Boltzmann constant 

domain (area or length) of a conduction element (see appendix A) 
3 a = 9/9x a , a = 1, 2 

9g = 9/9xg, £ = 1, 2 

Range of indices: 

I , J, K, L,N 1 to m (the number of nodes in the element) 

i,j,k,£,n Ito r (the number of reduced-system degrees of freedom) 

1,J,L 1 to n (the total number of temperature degrees of freedom) 

a,$ 1 to 2 

Superscripts: 

e element domain 

h convection 

k conduction 

o linear conduction 

r radiation 

1,2 normalized nonlinear conduction 


32 







1800 — 


1200 


X 

- Full system (190 D.O.F.) 
5 terms \ 

O 

6 terms 1 

Taylor 

A 

7 terms | 

series 

□ 

8 terms > 

1 


□ 

□ □ 


□ 


□ 



-600 L 

C 


D 


Figure 3.- Accuracy of solutions obtained by using Taylor series expansion for square plate with 

temperature-dependent thermal conductivities, q = 300 K. 
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Figure 4.- Accuracy of solutions obtained by using different sets of 
basis vectors for square plate with temperature-dependent thermal 
conductivities . 
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Figure 5.- Accuracy and convergence of solutions obtained with different sets of basis vectors 
(including the constant vector) for square plate with temperature-dependent thermal 
conductivities . 
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Figure 7.- Accuracy of solutions obtained by using different sets of basis vectors for cylinder having 
an eccentric hole and temperature-dependent thermal conductivities . 
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Figure 10.- Segment of the Shuttle wing bay 
used in the present study. 
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Figure 11.- Accuracy of solutions obtained by using single-parameter and three-parameter 
reduction methods for segment of the Shuttle wing bay. 
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Figure 12.- Accuracy and convergence of solutions obtained by single-parameter and three-parameter reduction 
methods (including the constant vector) for segment of the Shuttle wing bay. 
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